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ABSTRACT 

Numerical  solutions  are  obtained  for  two-dimensional,  unsteady 
potential  flow  about  a  flat  plate.   Two  types  of  flow  are  investigated: 
small  amplitude  oscillation  about  a  spanwise  axis  and  impulsive  accelera- 
tion of  a  nonoscillating  plate.   Pressure,  forces  and  moments,  and  kine- 
matics of  the  vortex  wake  are  calculated. 

In  this  work  vortex  sheet  shedding  from  the  trailing  edge  is 
included.   In  future  efforts,  vortex  sheet  shedding  from  trailing  and 
leading  edges  are  contemplated. 


NOMENCLATURE 

c  -  chord  length 

D '  -  drag  per  unit  span 

(i,  j,  k)  -  unit  vectors  in  (x- ,  y- ,  z-)  directions  respectively 

L'  -  lift  per  unit  span 

M  -  number  of  time  increments  taken 

N  -  number  of  evenly  distributed  bound  vortices 

A 

r;  -  position  vector 

A 

s  -   arbitrary  running  variable 

t  -  time 

T.V.S.  -  trailing  vortex  sheet 

A 

U  (t)  -  free  stream  velocity 

A 

U  -  final  free  stream  velocity 

A  t\ 

V(rQ)  -  field  velocity  at  rQ 

V.S.  -  vortex  sheet 

AAA 

(x,  y,  z)  -  absolute  reference  system 

A  A 

Y(s)  -  vortex  strength  distribution 

A 

T  -  free  vortex  strength 

9  -  amplitude  of  oscillation 

AAA 

(5  >  7\>    Q  -  body  fixed  reference  system 

CS?>  £7] »  £r)  "  unit  vectors  in  (§- ,  T|- ,  C-)  directions  respectively 

w  -  angular  frequency 


ii 


Nondimensional 

quantities 

A 

ci    =  L'/%p  U  * 

A     t 

cd   =  D'Ap  U«,| 

c 

c 

A 

r     =  r/c 

A 

s      =   s/c 

AAA 

(x,   y,    z)   =    (x,   y,    z)/c 

AAA 

(§,    Tl,    O    =    (?,    T],    Q/c 

A  A 

Uro(T)   =  U^CO/U 

A        A 

U)    =    UJC/Uoof 
A      A 

Y  =  Y/Uoo 

A       A 

r  =  r/uOTf  c 

A 

T    =    t    U00f/C 


Subscripts 

(x,  y)  -  component  in  (x- ,  y-)  direction 
(§»  Tl)  -  component  in  (§-,  T]-)  direction 
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I .   INTRODUCTION 

Present  generation  computers  make  numerical  approach  to  unsteady, 
two-dimensional  subsonic  flow  attractive.   In  this  study,  two  types  of 
flow  about  a  flat  plate  in  nonlinear  motion  are  investigated.   Although 
having  different  objectives,  related  studies  have  been  carried  out  by 
J.  P.  Giesing  [l],  N.  D.  Ham  [2~]y   and  F.  H.  Abernathy  and  R.  E.  Kronauer  [3] 

Observations  show  that  unsteady  flow  about  a  lifting  body  is 
accompanied  by  a  vortex  sheet  being  shed  from  the  trailing  edge.   In  the 
present  study  of  unsteady  flow  about  a  flat  plate,  the  plate  is  also  re- 
placed by  a  bound  vortex  sheet.   Both  the  bound  and  the  free  vortex  sheet 
are  then  replaced  by  a  finite  number  of  discrete  point  vortices;  thus, 
the  integral  equation  of  potential  theory  representing  the  field  velocity 
is  reduced  to  a  finite  difference  equation.   With  appropriate  application, 
this  finite  difference  equation  is  used  to  determine  the  positions  of  the 
shed  free  vortices  as  a  function  of  time  and  to  express  the  boundary 
condition  at  any  instant  of  time.   In  regard  to  the  boundary  condition, 
it  is  stipulated  that  the  strengths  of  the  vortices,  both  bound  and  shed, 
are  such  that  the  trailing  streamline  is  tangent  to  the  trailing  edge. 
This  satisfies  the  Kutta  condition.   Together  with  Kelvin's  theorem, 
it  is  sufficient  to  determine  uniquely  the  strengths  of  all  vortices  at 
any  given  instant  of  time.   Steady  flow  is  treated  as  a  degenerate  case 
of  the  above. 


The  kinematics  of  the  shed  vortex  sheet  and  strength  distri- 
bution of  the  airfoil  and  free  vortices  are  determined.   Calculation  of 
such  important  physical  parameters  as  stability  derivatives  follows 
immediately. 

II.   THE  BASIC  EQUATIONS 

A  potential  flow  model  is  used  for  unsteady  flow  about  an 
infinite  flat  plate.   Laplace's  equation  governs  the  velocity  field 
and  pressure  can  be  calculated  from  Bernoulli's  equation.   From  Kelvin's 
theorem  the  change  of  circulation  about  the  flat  plate  is  balanced  by 
the  shed  vortex  sheet (s)\.      Finally,  it  is  assumed  that  the  Kutta  condi- 
tion is  satisfied;  this  condition  is  sufficient  to  determine  the  rate 
of  vortex  sheet  shedding. 

Using  thin  airfoil  theory,  the  flat  plate  is  replaced  by  a 
vortex  sheet  where  vortex  strength  distribution  is  adjusted  to  make  the 
flat  plate  a  streamline  of  the  flow. 

Consider  an  arbitrary  flow  field  containing  a  vortex  sheet 

A  A 

with  vortex  strength  distribution  Y(s).  Vorticity  is  considered  positive 

A 

when  inducing  counterclockwise  circulation,   r  is  the  position  vector 


1.   Bluff  body  flow  is  stipulated  to  be  accompanied  by  vortex  sheets 
shed  from  the  leading  and  trailing  edges. 
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on  the  sheet  while  s  is  the  running  variable  along  the  sheet.   The  field 


velocity  at  r   is  then  given  by 


f 


A   A 


A   A 

Vfr..fc)  = 


T(s,t) 


(r    -  r) 


v.s. 


v.s. 


2fr 


k      x 


^ds    +   U(i)  L 

CD 


r    -  r 


Tfs,t) 


<; 


2* 


-i  (y. -y)   *  j(x.-x) 


(x.-x)    ^(y0-y)2Jj 


(i) 


"ds  -h        U   (t)  L 


The  same  equation  written  nondimensionally  is 


Vfr.,t)  = 


2fr 


< 


V.3. 


-i(y.-y)  +  jfx.-x) 


<xo-x)%(yc-y)2 


"ds    +   U  ft)  L 

(2) 


Hereafter  we  will  omit  specific  reference  to  the  temporal  argument. 

A  description  of  the  geometry  of  the  field  is  appropriate. 
The  absolute  reference  system  is  assigned  coordinates  (x,y,z)  with  their 
origin  at  the  geometrical  center  of  the  flat  plate.   In  the  present  study 
the  center  of  the  flat  plate  is  therefore  fixed  but  the  plate  is  free 
to  rotate  in  an  arbitrary  manner,   i,  j,  and  k   are  unit  vectors  in  the 
x,  y,  and  z  directions  respectively.  A  body  fixed  coordinate  system 
(§  >   'H*  O  with  unit  vectors  er,  e- ,  and  er,  also  has  its  origin  at  the 
geometrical  center  of  the  flat  plate. 
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Obviously,  any  point  in  the  two-dimensional  space  considered 
may  be  described  by  either  coordinate  system.  Suppose  £  is  the  vector 
from  the  origin  to  the  point  of  interest.   Then 


r 


I  X   + 


jy 


or 


Furthermore,  it  is  easily  seen 


Conversely, 


x  a     cos  e  §  -    sin  e  n 

y  =       sin  0  g    -f     cos  O  "^ 

£  =     cos  0X    +    sin  e  y 

q  »    -  sine  x   +    cose  y 


(3) 


(4) 


(5a, b) 


(6a, b) 


If  we  decompose  the  integral  into  the  plate  portion  of  bound  vorticity 
and  the  free  vortex  sheet  portion,  we  can  write  instead  of  (2) 


1/2  f  "1 

/  T(§)        -L(y.-sine£)  +  j(xo-cose|) 

/  —  i   "* — a    r  ^s 

*2  2^      (Xa-cose$)  +(yo-sineg)2  J 


/ 


T(S>  J-i(y.-y) 


T.V.S. 


+  ifx.-  x) 


25v     [    (x,-x)2  f  (y.-y)1  J 


(7) 


U   L 

Oft    .»-» 


6  - 


A  useful  result  obtained  from  (7)  is  the  velocity  normal  to 
the  plate  at  some  point  §   on  the  plate  (T]o  =  0—  ) .   Thus 


1/2 


Vn  (g    ,  0!) 


/0    \ 

-1/2 


?(£)      d£ 


S.-S 


/»  T(s)    T      cosG  (X-  X)  +   sine  (  ya  -  V) 

/    —     J      i s.  j 

*.     2k      [         (x0-x)2  +  (y0-y)2  J 


T.V.S. 
whereas  the  tangential  velocity  is 


-      U      Sin  G 

00 

(8) 


Vg  («..0J)   =    e^  •  V(|.,oU 


+       —  +     U    cos  © 


T(s)  Sine  (xo-  x)  -   cos  8  (y   -  y) 


2fr 


T.V.S. 


( *.  -  x  ) ;    *   <  y0  -  y ) 


ds 


(9) 


Now,  for  a  flat  plate  whose  angular  position  is  given  as  a  function  of 
time,  9(t),  the  nondimensional  boundary  condition  to  be  satisfied  is 


V,  (  U    =     £ 


de 


dfr 


--<  §   < -+ e 

2    °   2 


(10) 


9  (t)  is  an  arbitrary  input  to  the  calculation.  By  extending  the  boundary 
condition  an  indefinitely  small  distance  downstream  of  the  trailing  edge, 
the  Kutta  condition  is  satisfied;  i.e.,  the  trailing  streamline  is  tangent 
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to  the  trailing  edge.   To  determine  the  strength  of  shed  vorticity, 
one  makes  use  of  Kelvin's  theorem  where 


-  {  /T(s)  ds 


-   0 


(ii) 


The  limit  of  integration  can  include  any  portion  of  the  free  vortex 
sheet  in  a  Lagrangian  sense  or  it  can  include  the  airfoil  and  any 
contiguous  portion  of  the  vortex  sheet. 

III.   THE  FINITE  DIFFERENCE  SYSTEM  OF  EQUATIONS 

The  numerical  scheme  to  arrive  at  solutions  for  the  airfoil 
and  free  vortex  strength  distribution  and  the  kinematics  of  the  free 
vortex  sheet  will  now  be  discussed.   As  a  reference,  the  flow  chart  (see  Ap- 
pendix) will  be  instructive.   The  scheme  is  basically  two-part.   Knowing 
the  strengths  and  positions  of  all  vortices  at  t  =  t,  ,  the  velocity 
field  is  calculated  and  the  new  location  of  the  free  vortices  are  calcu- 
lated at  t  =  t,  +  At .   The  circulation  strength  of  each  free  vortex  is 
conserved.   Then  the  unknown  strengths  of  the  bound  vortices  and  most 
recently  shed  free  vortex  are  calculated  by  satisfying  boundary  conditions 
as  described  by  the  finite  difference  counterpart  of  (8)  and  (10). 

The  vortex  sheet  representing  the  airfoil  will  be  approximated 
by  a  finite  number,  N,  of  evenly  distributed  bound  vortices.   Bound 
vortex  positions  are  indicated  by  circles  in  Fig.  1  (where,  as  an  example, 
we  choose  N  =  5).  The  distance  between  the  vortices  is 


-1  (12) 

*      N-  1 
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The  matching  points,  5  . ,  at  which  the  boundary  conditions,  equation  (10), 
are  satisfied  will  be  halfway  between  consecutive  vortices,  so  that  as 
N  -»  co,  the  normal  velocity  will  approach  the  principal  value  of  the  integral 
given  by  (8).   Thus,  the  finite  number  of  evenly  distributed  bound  vortices 
will  be  placed  at 

J"1 


g.  =  -  0-5  + 
J  N-  1 


(j-  1.2 ,N) 


(13) 


while  the  matching  condition  points  will  be  located  at 


§.    =  -  0-5   + 


i  -  0-5 


U-i.2 .N-l) 


(14) 


N-l 

The  shed  vortex  sheet  will  be  approximated  as  discrete  doubly-infinite 
free  point  vortices  shed  from  the  trailing  edge.   A  new  vortex  is  shed 
at  each  increment  in  time  and  its  strength  is  uniquely  determined  in 
satisfying  the  boundary  conditions  (equation  10).   Once  the  strength  of 
each  vortex  is  determined,  it  will  remain  constant  for  all  time. 

To   start  the  calculation,  let  us  suppose  that  the  vortex 
sheet  distribution  is  known.  The  rate  of  change  of  position  of  the  i— 
free  vortex  during  a  time  increment  At  is  given  by  the  finite  difference 
counterpart  of  (7).  Replacing  the  vortex  sheet  strength  yLs   by  F.  at 


x,y  -*j.  7j 


V 


Xl 


_;    j 


y ■  -  sin  S  £. ) 


K 


r.i  2*      {  (xt-  case  £  )*+  (y.  -  sine  £.  ) 


(15a) 


1    r.    f 

>     -J  < 

2^ 
i-i 


(Vi-Vj) 


L'Xi-xo^fyi-yj)'] 


y  +    u 


oo 
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V  .    = 


N  r 

2     -J    <        — 

.  x  2*    L  (x.-  coseg.r  +•  (y-  sineg.)zJ 


f  X.  -    cosQ  £.) 
i 1JL_ 


Ag 


M    r 

2 J 


(x.-xj) 


1 


(15b) 


Then  the  new  position  of  each  free  vortex  may  be  determined 


X.    (t  ¥  At)   -    X.  (t)      +   V    At 
i+  1  I  xi 


y.   (  t  +at)  -   v.  (fc)  +-  V  •  At 

'4*1  7i  yi 


(16a, b) 


We  have  found  it  convenient  to  index  by  unity  each  free 
vortex.   Thus  the  newly  shed  free  vortex  will  always  be  identified  by 
i  =  1.   Having  established  the  positions  of  the  free  vortices  according 
to  (15)  for  t  =  t,  +  At,  combining  (8)  and  (10)  the  boundary  condition 
to  be  satisfied  at  that  instant  of  time  may  be  written  in  a  finite 
difference  sense  as 


N 


T. 


*£ 


i 

J      J      oi 


i    ii 


M 


-I 


f    ii 


(17) 


-  U  sin  e  -  § 


01 


d9 
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with  the  right-hand  side  completely  specified.   (Vf  . .)„  is  the  component 

t*Vi 

of  velocity  normal  to  the  surface  at  §  ,  induced  by  the  i —  free  vortex. 

J  oi 

It   is   given  by 


f  ii  J 


r 


£  .     -   (  cosS  x.   f  sine  y.) 
°i  J  'J 


2ft      Cg.  sine-  y-  )2  +  (|    cos©  -  x.)2 


(18) 


It  may  be  shown  that  Kelvin's  theorem  (11)  may  be  written  as 


H 


■  1  j  * 


i  ^  t.  ag 


1  j-i  J 


t    - 1  ♦'At  t, 

Equations  (17),  (18),  and  (19)  may  be  combined  and  written  as 


1    a..    b; 

■    1     J      J 

J-1 


c. 
1 


(t  -  1,2 ,NM) 


where 


(19) 


(20) 


M 


C 


2  KJ 

:.2         ftJ    ' 


-   U     sine    -    g. 


dO 
dt 


(21) 


N*l 


j"1 


t  •  k 


and, 


(22) 


(j-  i.2,....,N) 


(23) 
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H+l 


(24) 


The  coefficient  matrix  follows  readily 


i     n 


2*  *j  -  fi« 


i  -  1,2,....,  N 


I  J  -  J'2 N  J 


I    (25) 


i  N*l 


f  ii  J  1 


(  i  -  1,2 ,  N  )  (26) 


Finally  from  (18) 


N+-1 


a  s 


(j-  1,2, N) 


(27) 


a  =i.o 


(28) 


Numerical  representation  of  the  boundary  condition  and  Kelvin's  theorem 
is  now  complete.   Solution  is  obtained  by  inversion  of  the  coefficient 
matrix.    The  equations  were  programmed  in  Fortran  IV  for  the  IBM  360-67 
(see  Appendix)  . 
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IV.  CALCULATED  RESULTS 

The  solution  for  the  field  about  a  flat  plate  in  steady  uniform 
flow  is  first  considered.  A  comparison  is  made  in  Fig.  (2)  between  the 
vortex  strength  distribution  determined  by  the  finite  difference  technique 
and  that  predicted  by  the  exact  integral  (see,  for  example,  reference  [4]). 
Greater  agreement  is  achieved  by  increasing  N,  the  number  of  bound  vortices 
defining  the  flat  plate.   This  is,  as  noted  previously,  due  to  the  fact 
that  as  N  -»  °°,  the  solution  from  the  finite  difference  technique  approaches 
the  principal  value  of  (8). 

The  case  of  impulsive  acceleration  from  rest  was  considered 
with  8  =  0.1  rad  (5.73  ).   The  resultant  shed  vortex  sheet  is  depicted 
in  Fig.  (3)  for  T  =  2 .0  and  T  =  3.0.   The  vortex  strength  distribution  on 
the  plate  approaches  that  of  steady  state  asymptotically  with  time  while 
the  strength  of  successive  shed  vortices  becomes  negligible.   Fig.  (4) 
shows  ci  (normalized  on  angle  of  attack)  vs.  T.  As  expected,  it  approaches 
the  analytical  value  of  2tt  asymptotically. 

Harmonic  oscillation  in  pitch  was  then  considered.  9(t)  was 
defined  such  that 

e(T)  =  ©0  sin  oot 

The  important  physical  parameters  of  9Q  and  w   were  set  equal  to  0.1  and 
0.1,  1.0,  2.0,  and  3.0  respectively.   The  number  of  bound  vortices 
defining  the  plate,  N,  was  set  to  11  and  the  maximum  number  of  time  steps 
taken  was  150.   At  was  fixed  so  that  the  longitudinal  spacing  of  shed 
vortices  would  be  of  the  order  of  A§ . 
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Figs.  (5-10)  show  the  resultant  shed  vortex  configurations 
for  the  specified  frequencies  at  selected  elapsed  times.   Since  the  wake 
is  approximated  by  discrete  free  vortices,  curves  drawn  sequentially 
through  the  points  yield  an  approximation  to  the  shed  vortex  sheet. 
As  distance  downstream  increases,  definition  of  the  sheet  becomes  more 
difficult. 

The  short  period  frequency  of  a  stable  aircraft  at  cruise 
speed  corresponds  to  a  reduced  frequency  of  the  order  of  0.1.  As  seen 
in  Fig.  (5),  the  resultant  wake  has  negligible  vertical  deflection. 
This  reaffirms  the  assumption  that  the  wake  is  linear  for  most  practical 
purposes  in  unsteady  aerodynamics.   Comparison  of  Figs.  (6-10)  demonstrates 
the  increasing  nonlinearity  of  the  wake  with  increasing  frequency. 

An  interpretation  is  made  in  Fig.  (10)  of  the  vortex  sheet 
for  U)  ss  3.0  and  T  =  10.0.  At  a  sufficient  distance  downstream,  the  wake 
rolls  up  periodically  into  concentrations  of  free  vortices.   The  spacing 
ratio  (that  is,  the  ratio  of  the  vertical  displacement  of  the  "centers  of 
gravity"  of  cloisters  of  vortices  of  opposite  sign  to  the  horizontal 
distance  between  cloisters  of  like  sign)  is  typically  0.07.   The  classical 
ratio  predicted  by  von  Karraan  [5]  for  a  stable  configuration  of  double 
infinite  rows  of  point  vortices  of  alternating  sign  is  0.281. 

Coefficients  of  lift  and  drag  were  calculated  by  applying  the 
Kutta-Joukowsky  Law  locally  to  each  bound  vortex  defining  the  flat  plate. 
Figs.  (11,  12)  present  these  values  for  w  =  3.0  as  a  function  of  elapsed 
time  t.  A  typical  value  for  ci  at  T  =  0.50  and  9  =  -5.71  is  .28. 
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The  steady  state  value  for  this  particular  angle  of  attack  is  .61. 
Cx  and  c^  are  periodic  with  time  with  negligible  phase  shift  relative 
to  the  harmonic  oscillation. 

V.   CONCLUSION  AND  DISCUSSION 

The  method  developed  herein  yields  good  results  for  the  unsteady, 
two-dimensional  flow  about  a  flat  plate  executing  two  separate  nonlinear 
motions:  impulsive  acceleration  art!  small  amplitude  oscillation  about  a 
spanwise  axis.   With  the  flow  field  specified,  such  significant  quantities 
as  stability  derivatives  are  obtained  readily. 

In  future  efforts,  bluff  body  flow,  such  as  large  amplitude 
oscillation  or  tumbling,  with  vortex  sheet  shedding  from  the  trailing  and 
leading  edges,  is  contemplated.   Preliminary  investigations  were  made. 
Once  again  the  potential  flow  integral  equation  representing  the  field 
velocity  (1)  is  replaced  by  a  finite  difference  counterpart.   The  boundary 
condition,  (8)  and  (10),  is  extended  an  indefinitely  small  distance  upstream 
of  the  leading  edge  and  downstream  of  the  trailing  edge.   (15,  17,  18,  19,  20) 
are  modified  accordingly.   Kelvin's  theorem  determines  the  rate  of  vortex 
shedding  and  the  two  Kutta  conditions  establish  uniqueness.   Initial  attempts 
to  establish  the  flow  field  were  unsuccessful.   These  attempts  suggested 
that  further  refinement  is  needed  in  the  numerical  model.   Consideration 
must  be  given  vortices  in  close  proximity  (possibly  in  the  form  of  the 
criterion  proposed  by  Ham  [2]).   Furthermore,  some  modification  may  be 
required  for  the  boundary  condition  at  the  extremities  of  the  lifting  surface. 
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Figure  2.   Steady  State  Numerical  Solution  Compared  to  Classical 
Flat  Plate  Theory. 
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APPENDIX 

OSCILLATING  FLAT  PLATE 

COMPUTER  PROGRAM 


SPECIFY   INITIAL  CONDITIONS, 
CALCULATE   STEADY   STATE  FIELD 


SPECIFY  U    (t),    9(t) 


) 


CALCULATE  VELOCITY  FIELD 


'k+i  =  'k  +  At 


CALCULATE  POSITIONS  OF  FREE  VORTICES 


CALCULATE  STRENGTHS  OF  BOUND  VORTEX  SHEET 
AND  FREE  VORTEX  MOST  RECENTLY  SHED  BY 
SATISFYING  BOUNDARY  CONDITIONS 


C 


PRINT 


) 


PUNCH  DECK,  OSCILLATING  FLAT  PLATE, 

REAL*4NXM,NYM 

DIMENSION  E(llfiL),  Ft  11),  NG(ll),  NP(ll),  h(U),  U(ll) 

DIMENSION  XB(ll),  YB(li),  XF(20C),  YK200),  XM11),  YM11) 

DIMENSION  GAMMAF(200),  XiM(ll) 

DIMENSION  A(12,12),  6(12),  C(L2),  L(12)f  NV(12),  VNK(IC) 

DIMENSION  TGAM1200),  ZX(2G0),  ZY(2GO),  XlB(ll) 

PI  =  6.2831d 

N  IS  THE  NUMBER  OF  EVENLY  DISTRIBUTED  8LUND  VORTICES  CN  THE  LIFT- 
ING SURFACE.  AN  IS  THE  INCREMENTAL  DISTANCE  IN  THE  FINITE  DIFFE- 
RENCE SCHEME. 


T  IS  THE  ELAPSED 

INCREMENT. 


TIME,  TT  IS  THE  INITIAL  TIME  AND  DELT  IS  THE  TIME 


1HETAZ  IS  THE  AMPLITUDE  OF  OSCILLATION  IN  RADIANS  WHILE  FREU  IS 
THE  REDUCED  FREQUENCY  OF  OSCILLATION. 


13G 


N=ll 

AN  =  l.GVFLUAT(N-i) 

Nl  =  N-l 

N2  =  N+l 

THETAZ  =  0.1 

FKEQ  =  0.1 

CELT  =  C.l 

T  =  0.0 

TT  =  -1.0 

TZ  =  T  -  TT 

00  13  0  J  =  1,  N 

X18U)  =    -0.5  +  AN 

X1M( J )  =  -0,5  +  AN 

CONTINUE 


*  FLOAT(J-l) 

*  (0.5  +  FLOATU- 


1)  ) 


THETA  IS 
IN  THETA 


n.E  MGOIFIED 
WITH  RESPECT 


FREt  STREAM 
TO  TIME. 


ANGLE  WHILE  CIHETA  IS  THE  CHANGE 


140 


ISC 


YTE  *  YTE) 


ThETA  =  THETAZ  *  SINtFRtO  *  TZ) 

UcGRcE  -  THETA  *  57.29583 

ulHETA  =  THETAZ  *  FRlO  *  CGS(FREQ 

COSTH  =  CUS(  THETA  ) 

SINTH  =  SIM  ThETA) 

XTE  =  0.5  *  COSTH 

YTc  =  0.5  *  SINTH 

XXX  =  SCRTIXTE  *  XTE  + 

NXM  =  -YTE/XXX 

NYM  =  XTE/XXX 

Ou  14  0  J  =  1,N 

XL>(  J)  =  X1BU) 

Yb( J)  =  X16( J) 

XM(J)  =  XIM(J) 

YMIJ)  =  X1M(J) 

CONTINUE 

DO  150  I 

DO  15  0  J 

E(  I, J  )  = 

CONTINUE 


*  TZ) 


*  COSTH 

*  SINTH 

*  COSTH 

*  SINTH 


=  1 ,  N 
=  1,N 
i./(Pi 


(FLGAT( J-i )  -  0.5) ) 


c 
c 
c 
c 
c 


c 
c 
c 

c 
c 
c 


DO  35  0  I  =  lfN 

(•III  =  NXM 
j-30   CONTINUE 
IRUW  =  1 

CALL  M INV( E,N,D,NG,NP ) 
CALL  GMPRD(EfF,H,NfN,IROW) 
DO  60  I  -  L,  t\ 
811)  =  hi  I  ) 
60  CONTINUE 

\hc    INITIAL  BOUND  VORTEX  STRtNGTH   DISTRIBUTION  iS  SPECIFIED. 

CONTINUE 

lhE  UNSTtADY  FINITE  DIFFERENCE  SOLUTION  FOLLOWS. 


C 
C 
C 

C 

C 
C 

c 
c 


C 

c 
c 

c 

C 

c 


00  10C0  M=l, 150 
T  =  T  +  DELI 
TZ  =  T  -  TT 

ThE  FREE  VGRTICtS  AKL  ASSUMED  ShhD  FROM  THE  TRAILING  cDGE. 

X  F  (  1 )  =  X  B  (  N  ) 
YFl  1)  =  YB{N) 

ThE  VELOCITY  FlLLD  IS  DETERMINED. 

DO  BOO  J  =  i,  M 

VLLOCITY  INDUCED  BY  THE  I  TH  FREE  VOKTEX  ORIGINALLY  ShhD  FRGN  The 
TRAILING  EDGE  ON  ThE  JTH  FREE  VORTEX  ORIGINALLY  SHED  FKCM  1 nl 
TRAILING  EDGE. 

MXJ]  =  CO 
VYTT  =  0.0 
IF(M.EQ.l)  oO  TO  501 

GAMMAF(l)  IS  AS  YET  UNKNOWN. 

DO  50  0  I  =  2,M 

A  VORTEX  CANNOT  ACT  ON  ITSELF. 


IF( I.EQ.J )  00  TO  5  00 
X  =  XF{ J  )  -  XF( I  ) 
Y  =  YFU)  -  YF(  I  ) 
XX  =  X  *  X 
YY  =  Y  *  Y 

DENOM  -  PI  *  (XX  +  YY) 

VXTT  =  ((-Y  *  GAMMAF( I ) )/D£NOM)  +  VXTT 
VYTT  =  {(  X  *  GAMMAFt I ) ) /DENOM)  +  VYTT 
500   CONTINUE 
501  CONTINUE 


VELOCITY  INDUCED  BY  THt  ITH  BOUND  VORTEX  UN  THE  JTH  FRt, 
UKIGINALLY  SUED  FROM  THE  TRAILING  EDGE. 


VCKTEX 


VXBT  = 
VYBT  = 


CO 

0.0 


THE    BOUND  VURTEX 
THE  SINGULARITY. 


AT  XB(N)  IS  EXCLUDEO  TO  ARTIFICIALLY  SUPPRtSS 


600 


DO  60  0  I  =  1,N 

IFt I.EC.N.AND.J.EQ.l) 

X  =  XFI J)  -  XB( i  ) 

Y  =  YF( J)  -  YB(  I  ) 

XX  =  X  *  X 

YY  =  Y  *  Y 

DENGM  =  PI  *  (XX  + 

yAQ\    =     (  (-Y  *  (5(1) 

VYBT  =  (  (  X  *  B(  I  ) 

CONTINUE 


oO  TO  600 


YY) 

*  AN)/DfcNG.M) 

*  ANI/DENCM) 


+  VXBT 
+  VYBT 


THE  POSITIONS  OF  THE  FREE  VORTICES  ARE  DETERMINED  FOR 
T  =  T(K)  +  DELT. 


oOC 


660 


230 


ZX(  J) 
ZY(  J) 
CONTI 
DO  8t> 
XF(  J) 
YF(  J) 
CGNTI 
1  HETA 
Di&RS 
ALPHA 
DTHET 
COSTH 
SINTH 
XTE  - 
YTE  = 
XLE  = 
YLE  = 
XXX  = 
NXM  = 
NYM  = 
DO  2  3 
XB(  J) 
Y  B  (  J ) 
XM(  J) 
YM(  J) 
CONTI 


=  ( 
=  ( 

NUE 

0  J 
=  Z 
=  Z 

NUE 
=  T 

E  = 

A  = 
=  C 

=  s 

0.5 

0.5 
-XT 
-YT 
SCR 
-YT 
XIE 

0  J 
=  X 
=  X 
=  X 
=  X 

NUE 


VXT1  +  VXBT)  *  DELT  +  DELT  +  XF(J) 
VYTT  +  VYBT)  *  DELT  +  YF ( J ) 

=  L,M 

XI  J) 
Y(  J) 


HE 
TH 
(D 
TH 
OS 
IN 

if 

a* 

E 

E 

Tl 

E/ 

/X 


IB 
I  ft 


TAZ  *  SIN(FREQ  *  TZ) 

tTA  *  57#29583 

EGREt  ) 

ETAZ  *  FREO  *  CGS(FREQ 

( THETA) 

( THETA) 

COSTH 

SINTH 


*  TZ) 


XTE 

XXX 

XX 

1,N 

(J) 

(J) 

(J) 

(J) 


*  XTE  +  YTE  *  YTE) 


COSTH 
SINTH 
COSTH 
SINTH 


CALCULATING  THt  COMPONENTS  UF  THE  ROW  VECTOk,  Cll) 


EACH  ELEMENT  OF  THE 
CONDI T ION  POINT,  AN 
STREAM  VELOCITY  AND 
THOSE  MOST  RECENTLY 
TFt:  LIFTING  SORFACE. 


ROrt    VECTOR  IS,  AT  A  PAKTICGLAR  MATCHING 
EXPRESSION  OF  THE  NGRPAL  CCHPCNENTS  CF  FREE 
VELOCITY  INOOCED  BY  FREE  VCRT1CIES  (EXCcRT 
SHEO)  PLUS  THAT  COMPONENT  CUE  TO  ROTATION  OF 


TOG 
101 


Co  20  0  I  =  i,N 

CALCULATING  THE  INDUCED  VcLOCITY  AT  THE  I TH  MATCH  I  NO  CONDITION 
POINT  DoE  TO  THE  JTH  FREE  VORTEX  ORIGINALLY  SHEO  FROM  THE  TRAILING 
ECGE,  WITH  THE  EXCEPTION  OF  THAT  ONE  MOST  RECENTLY  SHEC. 

VX  =  0.0 

VY  =  CO 

IF  (M.EQ.l)  00  TO  101 

00  10  0  J  =  2,  M 

X  »  XM{  I )  -  XF(  J) 

Y  ■  YM(  I)  -  YF(  J ) 

XX  =  X  *  X 

YY  =  Y  *  Y 

DENGM  =  PI  *  (XX  +  YY) 

VX  =  ((-Y  *  GAMMAF( J) )/DENCM)  +  VX 

VY  =  ((  X  *  GAMMAFU)  )/DENOM)  +  VY 

CONTINUE 

CONTINUE 

VNM(I)  =  (VX  *  NXM  +  VY  *  NYM) 

Cll)  =  NXM  +  VNM(I)  -  IXlMli)  *  DTHETA) 

CONTINUE1 

Component  doe  to  kelvin's  THEORtM. 


b2G 


SUMGAM  =  0.0 

DO  o20  I  =  i,N 

SOMGAM  =  6(1)  *  AN  +  SUMGAM 

CONTINUE 

CIN2)  =  SOMGAM 


CONSTRUCTING  THE  COEFFICIENT  MATRIX. 


THOSE  ELEMENTS  DOE  TO  THE  VELOCITY  INDOCcD  BY  THE 
AT  THE  JTH  MATCHING  CONDITION  POINT. 


JTH  Euu.NO  VORTEX 


bu  210  I 

=  ltN 

DO  110  J 

»  if  N 

Al  I, J  )  = 

1./IP1 

110 

CONTINUE 

210 

CONTINUE 

*  (FLOATl J-I  )  -  0.5)  ) 


|"HOSE  ELEMENTS  DUE  TU  THE  EFFECT  OF  THE  FREE  VCKTEX  MOST  KECENTLY 
SHED  FROM  THE  TRAILING  EDGE  ON  THE  I  TH  MATCHING  CONDITION  POINT. 


530 


DO  53  G  I  =  L,N 

X  =    XM(  I  )  -  XF( 1) 

Y  =  Y  M  (  I  )  -  Y  F  (  i  ) 

XX  =  X  *  X 

YY  =  Y  *  Y 

DENOM  =  PI  *  (XX 

VX  =  -Y/DENOM 

VY  =    X/DENQM 

At  I,N  +  1  )  =  -(NXM 

CONTINUE 


+  YY) 


*  VX  +  NYM  *  VY) 


THGSt  ELEMENTS  DUE  TO  CONSIDERATION  OF  KELVIN'S  THEORY 


550 


CO  55  0  J  =1,N 
A(N2, J)  =  AN 
CONTINUE 

A(N2,N2)  =  1. 


Inversion  uf  ai i , j  ) 


CALL  Mii\V(  A,N2,DfL|NV) 


INVERSE! A( I, J))  *  C( I )  =  B(  J) 
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Ml'lM 
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GM 
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NU 
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2F 

,2 

•  I 

t 

I 

X 
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E 
1 
1  + 


DlAfC?B,N2iN2f IROW) 

=  B(N2) 
50)  GO  TO  999 


0) 
3.2 

3) 

10. 

6) 

10. 

1  ) 

10. 

2) 

COG 


T,  FREQ, 
,  FIG. 5) 
XLc, 


ALPHA 


YLE 
3) 

XTE,  YTE 
3) 

(XF{ J)  ,  YF( J) ,  J  =  1,M) 
5) 


) 


=  IfM 

GAMMAF ( I ) 
F(  I) 
F(  I  ) 

=  1 ,  M 

1  )  -  TGAM{ I  ) 

ZX  (  1  ) 

ZY(  I  ) 


thesW688  L< 
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